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LATERAL DENSITY VARIATIONS IN ELASTIC EARTH MODELS 


FROM AN EXTENDED MINIMUM ENERGY APPROACH 
B. V. Sanchez 
Abstract 

The gravitational field of the earth shows the existence of lateral density inhomogeneities in 
the mantle and crust but it does not yield a unique solution for them, that is, there ate an infi- 
nite number of possible density perturbation fields which will match the observed gravitational 
field. 


Kaula’s minimum energy approach has been extended to include the nonhydrostatic gravita- 
tional potential energy and the density perturbation field has been obtained to degree and order 
8. The depth profiles for the density perturbation show a stratification with density excesses and 
deficiencies alternating with depth. The results indicate that the addition of the gravitational 
potential energy in the minimization process does not change significantly the conclusions based 
on results for the minimum shear strain energy case, concerning the inability of the mantle to 
withstand the lateral loading elastically. 


PRECEDING PAGE BL ANK NOT FfLMEL’ 


CONTENTS 


Page 

Abstract iii 

1. Introduction 1 

2. Differential equations 3 

3. The elastic and gravitational potential energies 9 

4. Numerical results and conclusions II 

References 13 

ILLUSTRATIONS 

Figure Page 

4.1 Surfaces of equal density perturbation at the surface of the earth, 

units of gm/cm^. Minimum shear strain energy case 15 

4.2 Surfaces of equal density perturbation at the surface of the earth, 

units of gm/cmL Minimum gravitational potential energy case 16 

4.3 Surfaces of equal density perturbation at the surface of the earth, 

units of gm/cm^. Minimum shear plus gravitational energy case 17 

4.4 Surfaces of equal radial displacement at the surface of the earth, 

units of meters. Minimum shear strain energy case 18 

4.5 Surfaces of e({ual radial displacement at the surface of the earth, 

units of meters. Minimum gravitational potential energy case 19 

4.6 Surfaces of equal radial displacement at the surface of the earth, 

units of meters. Minimum shear plus gravitational energy case 20 

4.7 Surfaces of equal density perturbation on an equatorial cross-section 

of the earth, units of gm/cm^. Minimum shear strain energy case 21 

4.8 Surfaces of equal density perturbation on an equatorial cross-section 
of the earth, units of gm/cm^. Minimum gravitational potential 

energy case 22 

preceding page blank NOT 


V 


ILLUSTRATIONS (Continued) 


Figure Page 

4.9 Surfaces of equal dens!./ perturbation on an equatorial cross-section 
of the earth, units of gm/cm^. Minimum shear plus gravitational 
energy case 23 

TABLES 

Table Page 

4.1 Global energy magnitudes, units of lO^ergs 12 


LATERAL DENSITY VARIATIONS IN ELASTIC EARTH MODELS 
FROM AN EXTENDED MINIMUM ENERGY APPROACH 

1. Introduction 

The external gravitational potential of the earth can be expressed in spherical hannonics and 
its components determined by analysis of satellite data. Such an analysis indicates the presence 
of lateral variations of density in the earth interior. However, it is not possible to define uniquely 
what these lateral density variations are by inverting only the gravity data. The question of lateral 
density variations and the internal stresses they produce has implications with respect to the possi- 
ble existence of mantle convection within the earth. Kaula (1967), Lambeck (1976), Phillips and 
Lambeck (1980) among others have provided extensive and detailed discussions of these ideas. A 
possible approach to the problem consists in using static earth-tide theory in a modified form to 
include interior density anomalies expressed in spherical harmonics, the parameters which define 
the anomalies are determined by imposing the condition of minimum strain energy and the solu- 
tion is constrained to satisfy the external boundary conditions provided by the gravitational poten- 
tial; the surface topography is also included in the problem as a loading surface density distribu- 
tion. Tliis is the method followed by Kaula (1963) and Arkani-Hamed (1970, 1972). It provides 
an estimate of the minimum shear stresses induced by the internal density anomalies and conclu- 
sions can be drawn concerning the capability of the mantle to withstand such stresses elastically. 
The results indicate that the stresses are probably excessive. 

Lambeck (1976) has pointed out that the nonhydrostatic gravity field implies the existence 
of free potential energy within the earth and that there must be a capacity for rearrangements of 
mass toward distributions of lower energy (assuming no new energy generation and no kinematic 
constraints). A study by Rubincam (1979) indicates that the maximum nonhydrostatic gravita- 
tional potential energy in the mantle and crust is of the same order of magnitude as Kaula’s esti- 
mate for the minimum shear strain energy. 
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In the li^t of the remarks above, it is not without a measure of interest to carry out an in- 
vestigation similar to Kaula’s but with the inclusion of the gravitational potential energy in the 
minimization process. The object of the study is to determine the distribution of lateral density 
anomalies and to compare it to the results for the minimum strain energy case to ascertain if any 
significant changes occur. 
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2. Differential equations 


The equations of motion governing the vibration of an elutic body are given by Love (1944), 
these equations and the Poisson's equation for the gravitational potential can be expressed in spher- 
ical coordinates and with the appropriate set of boundary conditions they will admit an eigenvector 
solution. The static solution is obtained by setting the frequency equal to zero. Using the trans^ 
formation of Alterman et. al., (1959) these partial differential equations yield the following set of 
six ordinary differential equations. The dot denotes derivatives with respect to r. 

y'* IA)T+ IBl d* (2.1) 

T “ (yj. y2* yj. ys. ys> 

The elements of the matrix (Al are given by Smylie and Mansinha (1971) and Israel et. al., (1973). 
They arc given below for the sake of completeness. The symbols Pg* M and X denote the unpei^ 
turbed density and the elastic parameters, gg stands for the gravitational acceleration and n is the 


degree of 

the deformation. 


c 

= 1/(X + 2/i) 



= -2Xc/r 


^12 

= c 


^13 

= Xn (n + 1 ) c/r 


^21 

= 4p(3X + 2 m) c/r2 - 

4Po«o/f 

A22 

* ~4m c/r 


A23 

= -n (n + 1) [2p (3X 

+ 2m) c/r2 - pggg/rl 

^24 

= n (n 1 )/r 


^26 



Ajl 

= -l/r 


A33 

» 1/r 


^34 

* 1/M 


^41 

» -2m (3X + 2m) c/r2 

+ Po8o/r 


A42 ■ -Xc/r 

A4J ■ 4 n (n 1 ) #i(X ♦ n) c/r* - 2 #i/H 

A44 »- 3 /r 

A 4 S "-Po/f 

A 5 , -4*Gpo 

A55 ■ I 

A«j * -4»G n(n + I)po/r 

A« » n(n + l)/r2 
^ ■ ~2/r 

The unlisted elements are equal to zero. The yj are radial functions defined as follows 

“r = yt(r)S„m 

Ufi - y3(r)0/a<?)S„„, (2.3) 

* y3(r)(l/sin0)(d/d^)S„n, 

where u^, and denote displacements in the respective coordinate directions and Sn„| is a 
surface harmonic of degree n and order m. The stress-strain relations yield; 

“ (R 1 (V1.y2.y3. ( 2 - 4 ) 

where the 5's denote stresses and the matrix [R] is given by Arkani-Hamed (1972), 

Rj, = I2p <K3X + 2p)/rl 


R31 

- R 21 



* ®nm 


*^22 

* ^nm 


•^32 

“R22 


P23 

* (2p/r) [b^ldd^ - n(n + 

D^cl S„„, 

R33 

- (-2p/r) + 2n(n 

1) (X +p) c) S„n, 


= (2p/r sin 0) Id^/dd • 

- cot 6 a/ddl s„„ 
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R44 • (8/W) S„m 

R„ - (I/sin6)(d/d^)S„„ 

The unlisted elements are equal to zero. Also 

^ ■ y$(r)S„m 

(2.5) 

VtiO ■ Vs - ^'GpoVi 

where \p denotes the perturbation in the gravitational potential and G is the gravitational constant. 


The density perturbations arise from two different causes, the perturbation caused by the 
deformation due to the loading and a perturbation due to unspecified effects such as thermal and 
chemical processes. That is, 

p = Po(r) + Pd(r, 0,0) ( 2 . 6 ) 

where Pg(r) is a nominal radial distribution (from seismic data), p^j is the density perturbation 
due to deformation and is the added density due to unspecified processes. It is assumed that 
p^ can be expressed in spherical harmonics, 

* L Z ‘’nm IPanli , cos m0 + p/„^^ j sin m0) (2.7) 

n m 

where P^m an associated Legendre polynomial. Furthermore it is assumed that the coefficients 

PAnmi Pven as polynomials in the normalized radius, 
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PAnmi ‘ «.«) 

r» 

where the d,in,ij are constants to be determined and R denotes the outer radius of the earth. 


With these assumptions the forcing term in equation (2.1) is qiecified, the elements of the 


matrix (B) are given by 


®2k 

®6k 

B«k 


k * 

- 4irCHr/R)k-i 
0 fi - 1,3, 4, 5- k 


1 to 6 
■ 1 to 6 


(2.9) 
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The vector includes the constants to be determined 

‘^ml5* (2.10) 

Equation (2.1) is applicable to the earth mantle and crust, the equations for the liquid core are 
obtained by assuming hydrostatic equilibrium conditions, they are given by 
y, • 4tGpoys/so ♦ y* 

y* ■ ln(n + l)/r* - l6»Gpo/gorl y, - {*wCpJio + 2/r)yj 
and 

yi “ Vslio (2.11) 

yj ' 0 

y, = (4y5 + ryj)/n(r. + Dgo 

y4 » 0 

The boundary conditions appropriate for the solution of the equations are specified by the spher 
ical harmonic expansions defining the surface topography and the gravitational potential. Let the 

surface topography be given by 

• n 

T = S 51 Pnm (Tnmi cosm^ + T „„2 sinm0) (2.12) 

n*o m«o 

and the gravitational potential by 

“ n 

V * X) X] Pnm (Vnml sinm^) (2.13) 

n*o m “0 

The surface boundary conditions are then 

yi(R) ■ -8o (T„n,|) 

y4(R) • 0 

(2.14) 

yj(R) ■ 

y^fR) ♦ l(n D/Rl y,(R) • 4»GT„„,j 

The initial conditions for the s<^ution of equations (2.1 1) were approximated by an analytical 
sdution (.f the surface load problem for a model of the earth having a homogeneous mantle and 
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liquid core, Zharkov (1967). Values of yj and were obtained at r ■ R/IO. Equations (2.1 1) 
were then integrated numerically up to the core-mantle boundary using the earth model M) due 
io Landisman et. al., (1965) as given by Israel et. al.. (1973). This procedure was followed for 
every degree and order. The boundary conditions adopted at the core-mantle boundary are 
those developed by Israel et. al., (1973) and Crossley and Ckibbins (1975), that is, continuity of 
the variables y 4 and yj only. 


lart the solution of equation (2.1) be given by 

INI'J 


(2.15) 


where (M) and (N) are matrices and'?"is a boundary conditions vector. It follows tha‘. 

|M| = |A1 (Ml. (M(r„,ro)l = (I) 

IN] = I A) INI IB), IN(ro,ro)l » [O] 


(2.16) 


where (I) is the identity matrix. At some point r, 

■yir) * lM(r, ro)l‘c‘+ (N(r,ro)jT (2.17) 


The set of ordinary differential equations has been cast into a transition matrix form to facilitate 
the application of a least-squares minimization procedure. 


Assume the existence of a function H which can be written as 



C.I8) 


where (L| is a matrix of the appropriate dimensions. The function H is to minimized subject 
to the following equations of constraint. 


IE) 



(2.19) 


where |E| is a matrix andTis a vector of fixed boundary conditions. It can be shown (Arley 
and Buch, 1950) that the solution 
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t ' ■ 




( 2 . 20 ) 



(L1*T (El [(El (Ll-T (ElT]->r 


minimizes H in a least«squares sense. 




3. The elastic and gravitationat potential energies 

The expression for the total shear strain energy of the mantle is given by Arkani-Hamed 
(1970), 

■R 


^nmi " 


J 


[PIT 


(3.1) 


where denotes the core radius and the matrix [P] has the following elements 


P2t “ P 


12 

P 3 I “ Pl3 
P32 ' P23 


P,J = :(3X + 2m) 2 c2/3r^ 

P,j = -2(3\ + 2M)c2/3r 
Pj 3 = - n(n + I)P,i/2 
P 22 = 2c2/3 
P 23 = - n(n + 1 )P, 2/2 

P 33 = 2n^ (n + 1)7 (3\7 + 4 p 7 + 6 #iX) c^/r^ - n (n + l)/r7 
P 44 = n(n + i)/ 2 ii 7 

The unlisted elements are equal to zero. 

Making use of equation (2.15) it is possible to rewrite equation (3.1) as follows, 

- R 


■’nmi 


(3.2) 


“ffl /■ .V [IMI.INIIMPI d,Q 

The expression for the gravitational potential energy of the mantle is given by Rubincam (1979), 

(3.3) 


Unmi = {-16ff2G/(2n + 1)1 f p(r)rn+> f X"+2p(x)dxdr 

•'Rc -'Rc 


where p(r) denotes the radial component of the density perturbation which has two components: 
the one due to deformation (Pq) and the one due to unspecified causes (p^). Making use of 
equation ( 2 . 8 ), 

UA^. = [- 16ff2G/(2n + 1)1 "dT /* lr^"^r)-r(r) - Re^^^Tr)-? (RcOl dr T 

•/p 


(3.4) 


where 


“F(D * (gi. g2.g3.g4. 85* 8«) 

8 k = (r/Ry^-VCk + n + 2 ) k = 1 to 6 

^r) = (/i./2.^3.^4.^5»^6^^ 

4 * (r/R)*^*' k * 1 106 

The radial component of the density perturbation field due to deformation is given by 

Pd = PoAr - (3po/3r)yi (3.5) 

where A^ denotes the radial component of dilatation, 

^ = yi + (2/r)yi - [n(n + l)/rl yj (3.6) 

Inserting equation (3.6) into equation (3.5) and making use of the differential equation for yj 
yields 

Pd (3-7) 

where 

= (qi.qj.Qa.O, 0, 0) 

q, = -PolAi, + ( 2 /r) + {\lp^)bpJbT] 

Q 2 “ “PqAi2 

C.3 = - Pol A, 3 - n(n + l)/r] 


Making use of equations (3.7) and (2.15), equation (3.3) yields 

f [[M).(Nl]dxdr «|lj 


(3.8) 
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4. Numerical results and conclusions 


Equations (2.16) were solved numerically in order to evaluate the expressions for the shear 
strain and gravitational potential energies given by equations (3.2), (3.4) and (3.8), then equation 
(2.20) was used to obtain the best estimate for the vectors'? and d where the function K was 
representative of three separate cases: 

i) minimization of the shear strain energy, H 

ii) minimization of the gravitational potential energy, H = + U^n,j. 

iii) minimization of the sum of the two, H = S„,„j + 

The nominal earth model used is the M 3 model of Landisman et. al., (1965) as given by 
Israel et. al., (1973), the density and elastic parameters are given in table form and a cubic spline 
interpolation was used to obtain their values as functions of the radius. 

The spherical harmonic expansion used for the topography is the equivalent rock set given 
by Balmino et. al., (1973). The expansion used for the gravitational potential is the Goddard 
Earth Model 6 by Lerch et. al., (1974) with the hydrostatic components of the second and fourth 
degree zonals subtracted out (Munk and MacDonald, 1960). 

The density perturbation field due to unspecified causes was obtained for each degree and 
order by means of equations (2.7) and (2.8), equation (2.7) yields the total sum which in this 
case was truncated at degree and order 8 . The surfaces of equal density perturbation evaluated 
at the surface of the earth are shown in figures (4.1), (4.2) and (4.3) for each of the three cases 
mentioned above. Figures (4.4), (4.5) and (4.6) present the surfaces of equal radial displacements 
at the surface. Figures (4.7), (4.8) and (4.9) show the surfaces of equal density perturbation 
evaluated on an equatorial cross-section of the earth. The magnitudes of the density perturba- 
tions and displacements at the surface is larger by about a factor of two for the case of minimum 
shear plus gravitational as compared to the minimum shear case. The same result is observed for 
the density anomalies as a function of depth. The minimum gravitational potential energy case 


yields much smaller density anomalies but unrealistically large displacements. The depth profiles 
for the density perturbation show a stratiFication with density excesses and deficiencies alternating 
with depth, this internal compensation has been mentioned by Phillips and Lambeck (1980) and 
others; on the other hand Arkani-Hamed (1970) refers to an oscillation along the radius as a nu- 
metical instability which he eliminates by introducing a weighting factor for the amplitudes. The 
minimum gravitational energy case yields abrupt lateral variations also. 

Table (4.1) gives the global energy magnitudes for each case. The minimum total energy is 
obtained when minimizing only the shear strain contribution. 

The results indicate that the addition of the gravitational potential energy in the minimiza- 
tion process does not significantly alter the lateral density perturbation field, as compared to the 
results obtained for the minimum shear strain energy case: certainly not in a way which would 
change the existing conclusions concerning the inability of the mantle to withstand the resulting 
stresses elastically. 


Table 4.1 

Global energy magnitudes, units of 10^ ergs 


DFGRFE 

Min. Shear 

Min. Grav. 


Grav. 

Shear 

Grav. 

Shear 

mam 

- 7.912 

0.178 

-0.086 

94.639 

z 

II 

- 5.966 

0.167 

-0.120 

204.381 

N = 4 

- 6.020 

0.566 

-0.237 

586.165 

II 

z 

-13.935 

0.225 

-0.321 

1062.685 

N = 6 

- 1.412 

0.070 

-0.174 

719,578 

N = 7 

- 0.672 

0.043 

-0.274 

1249.120 

N = 8 

Bsa 

0.029 

-0.324 

1307.191 

£ 

N 

-36.351 

1.281 

-1..539 

5223.762 

Total 

-35.070 

5222. 



Min. (Shear + Grav.) 


Grav. 


-27..S99 


-13.683 


-11.624 


- 7.443 


- 1.645 


- 3.510 


.623 


-68.130 


-55.108 


Shear 


5.640 


1.120 


1.039 


0.688 


0.062 


2.204 


.266 


13.021 































































References 



I 

I 

r 

V 

f 




I 


Alterman, Z., H. Jtrosch, and C. L Pekeris, I9S9, “Oscillations of the earth”, Proceedings of the 
Royal Society, A, 252. 

ArkannHamed, J., 1970, “Lateral variations of density in the mantle", Geophys. J. R. astr. soc., 

20 . 


Arkani-Hamed, J., 1972, “Density and stress distribution in the moon”. The Moon, 7. 

Arley, N., and K. R. Buch, 1950, Introduction to the theory of probability and statistics , John 
Wiley and Sons. 

Balmino, G., K. Lambeck, and W. M. Kaula, 1973, “A spherical harmonic analysis of the earth's 
topography”. Journal of Geophysical Research, 78, no. 2. 

Crossley, D. J.. and D. Gubbins, 1975, “Static deformation of the earth’s liquid core”. Geophys- 
ical Research Letters, 2. 

Israel, M., A. Ben-Menahem. and S. J. Singh, 1973, “Residual deformation of real e^rth models 
with application to the chandler wobble”, Geophys. J. R. astr. soc., 32. 

Kaula, W. M., 1963, “Elastic models of the mantle corresponding to variations in the external 
gravity field”. Journal of Geophysical Research, 68, no. 17. 

Kaula, W. M., 1 967, “Geophysical implications of satellite determinations of the earth's gravita- 
tional field”. Space Science Reviews. 7. 

Lambeck, K., 1976, “Lateral density anomalies in the upper mantle”. Journal of Geophysical 
Research, 81. no. 35. 

Landisman, M., Y. Sato, and T. Nafe, 1965, “Free vibrations of the earth and the properties of 
its deep interior”, Geophys. J. R. astr. soc., 9. 


13 


Lerch. F. J„ C. A. Wagner, J. A. Richardson, and J. E. Brownd, 1974, “Goddard earth models 
(5 and 6)", X-921-74-I45, Goddard Space Right Center. 

Love, A. E. H„ 1 944, A treatise on the mathematical theory of elasticity , Dover Publications. 

Munk, W. H., and G. J. F. MacDonald, I960, “Continentality and the gravitational field of the 
earth”. Journal of Geophysical Research. 65. no. 7. 

Phillips, R. J., and K. Lambeck. 1980. “Gravity fields of the terrestrial planets: long-wavelength 
anomalies and tectonics". Reviews of Geophysics and Space Physics. 18. no. 1. 

Rubincam, D. P.. 1979, "Gravitational potential energy of the earth: a spherical hannonic 
approach". Journal of Cieophysical Research, 84, no. Bll. 

Smylic, D. F., and L. Mansinha. 1971, “The elasticity theory of dislocations in real earth models 
and changes in the rotation of the earth". Geophys. J. R. astr. soc., 2 ^. 

Zharkov. V. N.. 1967. "Flementary theoiy of spheroidal vibrations of the earth". Academy of 

Sciences of the U.S.S.R. (Physics of the solid earth). 


14 


ircALriY 











Surfaces of equal density perturbation at the surface of the earth, units of gm/cm^ . Minimum gravitational potential 







of gm/cm^ . Minimum shear jrfus grawiUtional 









diq>lacefnent at the surface of the earth, units of meters. 






19 


Figure 4.5. Surfaces of equal radial displacement at the surface of the earth, uruts of meters. Minimum gravitational potential 
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